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TECHNICAL PUBLICATION 


BASELINE COMPUTATIONAL FLUID DYNAMICS METHODOLOGY FOR 
LONGITUDINAL-MODE LIQUID-PROPELLANT ROCKET COMBUSTION INSTABILITY 

1. INTRODUCTION 


Liquid-propellant rocket combustion instability research has been pursued with varying degrees of 
emphasis for more than 50 yr. These efforts have led to the successful development of many reliable 
engines, which have been mainly the result of much trial and error development testing based on a great 
deal of intuition and experience, and engine design remains today as much engineering art as science. A 
detailed understanding of the mechanism by which combustion instability occurs in these devices simply 
does not exist, and the occurrence of resonant combustion instabilities continues to be a major risk in the 
development of any new liquid-propellant rocket engine. 

Previous predictive methods have been based primarily on approximate analytical models, but 
computational fluid dynamics (CFD) capability has advanced to the stage where it can now be effectively 
used as a research and development tool. 1 Clearly, CFD methods have the potential to resolve the funda- 
mental chamber processes in great physical detail, yet these results are highly dependent on the validity 
and fidelity of the various physical submodels. For this reason, CFD, at this stage, is more valuable as a 
computational test-bed than as a practical design tool. 

Utilization of CFD methods for the combustion instability problem also requires careful consider- 
ation of numerical accuracy. For example, both the dissipative and dispersive characteristics of a particular 
scheme can have an important influence on numerical accuracy in unsteady flow problems, as demon- 
strated by Hsich.- Furthermore, stability limits can be greatly affected by the transmission and reflection 
of flow disturbance incident on the boundaries, and implementation of appropriate boundary conditions is 
not trivial. 

This Technical Publication develops and evaluates a computational method for the analysis 
of longitudinal-mode liquid-propellant rocket combustion instability based on the unsteady, quasi- 
one-dimensional Euler equations. The unsteady Euler equations in inhomogeneous form retain full 
hyperbolicity and are integrated implicitly in time using a second-order, high-resolution, characteristic- 
based, flux-differencing spatial discretization with Roe averaging of the Jacobian matrix. Combustion 
process source terms were introduced through the incorporation of a two-zone, linearized representation: 
(1) A two-parameter collapsed combustion zone at the injector face, or (2) a two-parameter distributed 
combustion zone in which interphase transport is derived from a Lagrangian treatment of the propellant 
spray. The method is evaluated against a simplified analytical solution based on linearized small distur- 
bance theory, and the numerical methodology is then exercised on a generic combustor configuration 


using both collapsed and distributed combustion zone models with a short-nozzle admittance approxima- 
tion for the outflow boundary. 

The resulting baseline open-source CFD code is able to serve educational/pedagogical needs in its 
current form and could eventually be developed into a practical research and development tool by extend- 
ing the methodology to multiple dimensions, incorporating realistic physical submodels, and engaging in 
a painstaking validation effort. Ultimately, CFD methods may be of most practical utility as computational 
test-beds for investigating and studying the underlying physico-chemical mechanisms associated with 
liquid-propellant rocket combustion instability. 
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2. THEORETICAL AND MATHEMATICAL FORMULATION 


2.1 Governing Equations 


The gas phase is treated using an Eulerian description with appropriate source terms for interphase 
transport coupling with the spray. Thus, the problem reduces to solving the unsteady, quasi-one- 
dimensional Euler equations, which may be written in conservative form using matrix notation as 


dt dx 

where the conserved variable vector, Q, and the convective flux vector, E, are defined by 


( 1 ) 


and 


pS 


Q = 


puS 


peS 


( 2 ) 


E = 


puS 

| pu 2 + pjs 

( pue + up)S 


( 3 ) 


The source vector, El, contains contributions due to the variable cross-section area and interphase 
transport source terms. This vector is denoted as the sum of two components: 


H =H 1 + H 2 . 


( 4 ) 


For single-phase flow, only H j need be considered to account for the variable area effect on the momen- 
tum equation. It consists of a single component: 


H, 


0 

dS 

i — 

dx 

0 


( 5 ) 


For two-phase flow, the source vector H 0 must be included to account for interphase transport effects on 
the gas-phase conservation equation. It has the form 
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(6) 


H, 


Sw 

SwUp + D 


Sw 




+ Up + SuD 


Here, vv represents the propellant vaporization/combustion rate per unit volume. The perfect gas equation 
of state p = pRT = ( 7 - 1 ) e- pu 2 / 2 completes the governing system, where e = C v T + irll repre- 
sents the total internal energy. 


This inhomogeneous form of the Euler equations retains full hyperbolicity. Thus, the Jacobian 
coefficient matrix, A = dE/dQ , has real eigenvalues (u, u + c, and u- c with the speed of sound defined 
by c = Jyp/p). Following conventional practice, A may be diagonalized such that 

T~ l M~ l AMT = A , (7) 

where M = dQ/dV is a transformation matrix from the conservation variables Q to the primitive variables 
V, T is a similarity transformation matrix formed from the eigenvectors of the primitive variable Jacobian 
coefficient matrix M _1 AM, and A = diag {u, u + c, u- c} is the diagonal eigenvalue matrix. 

For convenience, positive and negative eigenvalue matrices may be defined as A ± , where A + con- 
tains only positive eigenvalues and A~ contains only negative eigenvalues. It therefore follows that 

A = Mt(a + +A~\t~ 1 M~ 1 = A + +A“ and \A\ = MT (a + - A~)t~ 1 M~ 1 = A + - A~ . 

2.2 Numerical Scheme 


Integration of the Euler equations is accomplished through temporal and spatial finite- 
difference discretization with appropriate linearization. For spatial discretization, a conservative formula- 
tion was desired that could yield good spatial resolution while avoiding spurious numerical oscillations. 
In selecting a time-integration method, the central concerns were obtaining adequate time-accuracy and 
good numerical stability. These attributes are central to the description of unsteady flows typically encoun- 
tered in a combustion stability analysis. Therefore, both the dissipative and dispersive characteristics of a 
numerical scheme are important points of consideration. 

An assessment of numerical techniques for unsteady flow calculations has been carried out by 
Hsieh. 2 Based on his survey of various spatial discretization schemes, sixth-order central differencing of 
the convective fluxes and fourth-order artificial dissipation yielded very good all around performance. He 
also explored characteristic-based schemes and found that second-order upwind scheme (2UP) flux dif- 
ferencing using Roe-averaged Jacobian coefficient matrices and nonlinear flux limiters to achieve total- 
variation-diminishing (TVD) conditions performed very well with regard to shock capturing but entailed 
higher dissipative errors in multidimensional problems. This characteristic based spatial discretization 
scheme was adopted, along with a generalized implicit time integration method, for the present develop- 
ment. Although not quite as effective as an explicit Runge-Kutta multistage approach, implicit time inte- 
gration provides adequate resolution and stability and is very robust. 
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2.2.1 Discretization 


Development of the numerical scheme begins with a generalized implicit formulation for the Euler 
equations in the form 


AQ l ] + e 


At 

Ax 




-i«+l 




= 0 , 


( 8 ) 


where 9 = 1 yields the Euler implicit (El) scheme and 9 = 1/2 yields the Crank-Nicholson (CN) (trape- 
zoidal) scheme. The first term is defined by AQf = Q'j + 1 - Qj, and the numerical flux, E, is detailed to 
first-order accuracy by 


p(D - - 
,/+l/2 - 9 


E j + E j + 1 


^1 /+ 1/2 ( Q J +1 Q j) 


(9) 


where A is the Jacobian coefficient matrix evaluated using the Roe-averaged quantities 3 ’ 4 


Pj+ 1/2 - Jpj+lPj 


( 10 ) 


1 j+l /2 



pj + “J+W 

'p/+l 

A 


'pj + 1 


( 11 ) 


h j+l /2 


h j\l 

pj + Vw 

'p./+i 

A 

fij + v 

'pj+i 


( 12 ) 


and 


'J+l/2 


(r-i) 


2 

, U J+l/2 

^ J+l/2 


(13) 


The delta form of the generalized implicit scheme is obtained by introducing AE n = E n+ ^ - E n 
and AH n = H n+l - H n = dH/dQAQ n : 


AQ'j + 9 


— ( AE n : 


Ax 


- AEt,n\-At ^ 

SQ 


j+l /2 ~ ^ j- 1/2 


AQ'j 


-RQ'I 


(14) 


where 


5 


(15) 


AE]+ 1/2 


2L 


AE j + AE j + 1 _ 1^1 /+1/2 ( AQ J+ l ~ AQ J ) 


and R\Qj is a residual explicit operator defined by 


R 



~ E i- 1/2 




(16) 


The variable cross-sectional area effect is included in the source vector Jacobian, but all interphase 
transport appears in the residual explicit operator. To arrive at the final form used for computation, substi- 
tute AEj = AjAQj and collect l ik e terms to yield 


- e —( A j-i + \ A \"- ll 2 ]\ A Qj-i 


2 Ax \ 


+ < 


SQ 


+e — |ur 


3 


2AxY 'J +l/2 


+ A 


if— 1/2 


AQ) 


+ id 


At 


-(a? 


2 Ax \ 


7+1 


+ A 


7+1/2 


> A Q n j+ i • 


(17) 


This first-order accurate scheme exhibits robust nonoscillatory behavior, but spatial resolution is 
poor. Direct substitution of second-order difference formulas does not improve matters either because this 
leads to the generation of spurious oscillations. To circumvent this difficulty, Harten has introduced the 
TVD property and derived sufficient conditions for constructing schemes that satisfy this property. 5 Such 
constructions have been carried out for both explicit and implicit time integration schemes. 6 ’ 7 The mecha- 
nisms currently in use for TVD conditions are based on some kind of gradient limiting procedure. We 
adapt the procedure based on a nonlinear flux limiter. The second-order convective flux vector is obtained 
by adding corrective terms to the first-order flux vector. A generalized form which may be specialized for 
second-order central difference (2CD), 2UP, or third-order biased upwind scheme (3UP) is given by 2 


p( 2) _p(l) , l ~ K 

h j+\l2~ h j+\l2 + 4 


AEj-H2 ~ AEj +3/2 


1 1 + K 

+ 


AEj +\/2 - AE i 


.7+1/2 


(18) 


where 


AE]+ m =(MT). m A& 


_/+1/2 7+1/2 


AE~j+]/2 - ( MT )j+ m Adj+i/2 


Ad j+y 2 = min mod 


AO j+ll2’PAOj_y 2 


(19) 

( 20 ) 
( 21 ) 
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( 22 ) 


Ao y+i /2 = min mod 


Ac 


,/+l/2^ zicr 7— 3/2 


4^1/2 = ^ -HW ( r_1 (2/+1 -Cy) • 
K + ,/2 = f f* + M); +1/2 ) y+1/2 (fij+1 - Qj) • 

and the minmod operator is defined by 

min mod [a;, y] = sgn(x)max[0,min{|x|,ysgn(x)}] . 


(23) 

(24) 


(25) 


The values of k corresponding to the alternative spatial discretization schemes are K = 1 for 2CD, 
K = — 1 for 2UP, and K = 1/3 for 3 UP. 

The constant, /3, is a compression parameter that is restricted to the range 1 < f5 <(3 -k)/(1-k) with 
1 3-6 when k = 1 . To maintain an efficient algorithm, some time accuracy is sacrificed by including the 
second-order numerical flux in the residual explicit operator only. In this way, the tridiagonal block struc- 
ture of the first-order flux is maintained; thereby, avoiding a pentadiagonal block system. 


2.2.2 Boundary Conditions 


Boundary conditions for the computational analysis are implemented according to the propagation 
of information along the flow characteristics. On a locally one-dimensional basis, for instance, each 
eigenvalue may be associated with a particular characteristic. If the eigenvalue is positive, the correspond- 
ing characteristic is right running. Conversely, if the eigenvalue is negative, the characteristic is left run- 
ning. When a characteristic runs out of the computational domain, the boundary condition depends on the 
internal flow field, and a numerical boundary condition must be formulated that depends on the computed 
solution. The number of unknown flow parameters minus the number of numerical boundary conditions 
gives the number of physical boundary conditions that must be imposed. The numerical and physical 
boundary conditions must be compatible with the system of equations for the characteristic variables. 

The numerical boundary conditions are formulated for an implicit scheme from the characteristic 
form of the equations with a selection matrix for choosing the outgoing characteristics. Also, one-sided 
differences are used to couple a boundary node to the interior. The resulting forms are as follows: 

Inflow boundary ( j = l) : 


l(t l M 1 
Outflow boundary (_/ = /): 


AQU +0 — AE" +m -At — 
■ J Ax/2 J+l12 dQ 


AQ, 


= -l[t l M ^rIq]) • (26) 
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AQ'j +6 


— AE"_ 1/2 
Ax! 2 J 1/2 


z\/ 


c)H 

dQ 


AQ, 


=- L ( 


T l M 1 




(fi 


(27) 


Here, 7'- 1 A/- 1 is introduced to obtain the characteristic variables and L is the characteristic selection matrix. 
For example, to choose the left-running characteristic only for subsonic inflow, L = diag {0,0, l}. In the 
case of supersonic outflow, there are three outrunning characteristics such that L - diag (l, 1, l}. If the 
outflow boundary is subsonic, on the other hand, there are only two right-running characteristics, and 
L - diag {1,1, 0j. By using one-sided differences, each numerical boundary relation has a two-node sup- 
port involving a boundary node and an adjacent interior node. 

* 

The imposed physical boundary conditions are imposed by letting B represent the specified 
boundary flow parameter. Then, Newton iteration implies B /?+1 = B" +dB/dQAQ n , and by enforcing 
B n+l -> B* , we obtain the constraining relationship 



AQ n j = B* - B n . 


(28) 


Both the linearized numerical and physical boundary relations are solved with the system of linear alge- 
braic equations for the interior nodes. 


2.3 Combustion Model 

A generalized schematic of a conventional liquid propellant rocket engine is shown in figure 1. 
Propellant combustion in the rocket chamber is modeled as a two-zone process consisting of a collapsed 
combustion zone at the injector face boundary and a distributed combustion zone extending from the 
injector face to the point where propellant reaction is completed. Thus, the inflow boundary for the com- 
putation is a transpiring reactive surface in which a small fraction of the propellant enters as a hot subsonic 
gaseous flow while the bulk of the propellant enters in the form of an unreacted liquid spray. Thus, the 
subsonic inflow boundary has only one outrunning characteristic, and two physical boundary conditions 
must be specified. 

The propellant combustion rate can be sensitive to pressure and velocity fluctuations, which per- 
mits the development of a feedback loop for combustion-driven instabilities. Here, simplified two-param- 
eter linear formulations are introduced to account for these sensitivities in both the collapsed and distributed 
combustion models. Furthermore, a two-parameter linearized spray atomization model has also been 
incorporated to reflect similar sensitivities associated with this process. The computational domain and 
the various submodel regions are illustrated in figure 2. 

2.3.1 Collapsed Combustion Zone 

Imposing a collapsed combustion zone at the injector face provides a mechanism for burned gas- 
eous propellant to enter the computational domain at the inflow boundary. Thus, from a physical perspec- 
tive, the inflow boundary may be viewed as a transpiring reactive boundary layer. Here, the characteristics 
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of this transpiring reactive boundary may be defined by the mean burned gas Mach number, M^, and the 
mean burned gas temperature, 7)-,. For convenience, we express the injector face combustion source as the 
sum of a steady mean flow quantity and an unsteady perturbation: 


m s = m s + m s , 


(29) 


where the mean mass transpiration rate corresponds to the steady solution for the particular rocket engine 
under consideration. 

Sensitivity of the transpiring combustion rate to pressure fluctuations and the time rate of change 
in these pressure fluctuations is accounted for through a simple phenomenological two-parameter linear- 
ized model in the form 8 


where the overbars symbolize mean quantities and primes represent perturbations. The second term on the 
right-hand side introduces a characteristic time scale defined as the ratio of the combustor length, l c , to the 
mean acoustic velocity, c. 

2.3.2 Atomization Zone 

The remainder of the propellant not consumed in the collapsed combustion zone is assumed 
to enter the computational domain in the form of liquid jets that completely atomize over a predefined 
length, I a . The resulting spray is characterized by a Sauter mean diameter (SMD), which is then used to 
compute the effective interphase transport terms. Because the atomization process is sensitive to both 
pressure and velocity fluctuations, the SMD can actually display a time-dependent variation that is directly 
coupled to chamber oscillations. To capture the underlying fundamental coupling process, a linearized 
two-parameter atomization model is adapted using the form 


where a a is the pressure sensitive index and (3 a is the velocity sensitive index for atomization. Numerical 
values for these parameters have been derived for various atomization models and are tabulated by 
Grenda et al. 9 

2.3.3 Distributed Combustion Zone 

The atomized propellant spray is represented as a collection of discrete computational parcels 
which are injected into the combustor at specified time intervals. Each parcel is identified with a group of 
physical droplets, the number of which is determined from the propellant injection rate and the mean drop- 
let size. The spatial distribution of the physical particles within each computational parcel is defined by a 



(30) 



u 


u 


(31) 
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probability density function (PDF) in space. The behavior of the physical droplets is determined by the 
history of their representative parcel, whose position corresponds to the mean of the PDF. The parcels are 
tracked through the combustor using a Lagrangian procedure, and knowledge of the shape and statistical 
parameters of the PDF permits calculation of interphase transport terms. 

PDF representation is utilized in order to minimize the numerical shot noise that can result when 
each parcel is described by a delta function point source and the ratio of parcels to physical droplets in the 
spray is too small. The numerical noise arising from a point source distribution function can be minimized 
by using a large number of parcels; however, such an approach can be extremely costly in terms of com- 
putational efficiency. 

Flere, we employ a rectangular uniform PDF shape that is fully defined given the mean position 
and half- width. The streamwise width of the PDF corresponds to the product of the injection velocity and 
a specified injection time interval, which is generally set to obtain a PDF width on the order of the com- 
putational cell width. In addition, the PDF shape and width remain fixed throughout a parcel’s lifetime. 


The mean position of each PDF is determined from the Lagrangian tracking of a single computa- 
tional parcel characterizing the behavior of the physical droplets in the group. Thus, parcel motion is 
governed by the the equations of motion, 


diip 

dt 


3 Pg C D 

4 Pp dp 


u - Up 


( 


U - Up 


) 


(32) 


and 


dxp 

dt 


(33) 


For the sake of simplicity, these equations are linearized by holding {^p g / Pp^yCj^ j d p j | 
over the computational time step and integrating analytically. This gives the updated parcel velocity as 


u-u, 


constant 


Up ~ Up Q exp 


-At 

T 


+ u\ 1 - exp 


-At 


(34) 


and the parcel position is obtained using explicit Euler integration 


% p ^ pfi ^ p,0^ • (35) 

The unsteady combustion rate for each parcel is assumed to be governed by the vaporization rate, 
which can be linearized about some mean steady state to obtain a vaporization-controlled combustion 
response function. For convenience, we express the local propellant combustion rate as the sum of a 
steady mean flow quantity and an unsteady perturbation: 

trip = dip + dip , (36) 
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where the instantaneous vaporization/combustion rate has the general form 

m p = m p ( K p + K u ) (37) 

and Kp and K u represent the pressure and velocity sensitivities, respectively. The mean mass burning rate 
of a parcel is governed by the mean burning rate constant, (5 , such that 


— 71 — 

nip — — ppdp/5 . (38) 

Note that the evaluation of the instantaneous propellant vaporization/combustion rate requires 
knowledge of the droplet temperature variation with time. The simplest feasible model for internal energy 
transport that can account for local gas-dynamic fluctuations is based on the rapid-mixing limit assump- 
tion. Here, the droplet temperature is assumed to be uniform due to strong internal circulation, but its time 
rate of change is allowed to vary according to the balance between the energy leaving the surface due to 
vaporization and the energy arriving due to heat transfer. Linearization of this model leads to a parcel 
combustion response function of the form: 9 


tn„ * p' 

-=- = u v — 

m p p 



(39) 


where the pressure- and velocity-sensitive indices, a v and p v „ are complex coefficients due to the phase 
difference that can exist between the vaporization and gas oscillations as a result of temporal oscillations 
inside the droplet. This phase difference can also be accounted for through the introduction of appropriate 
time lags when evaluating the fluctuating gas-dynamic properties such that the response indices are real 
numbers 


m p p u 

These combustion response time lags can have important effects on combustor stability and, in general, 
may not be neglected. 



2.4 Short-Nozzle Approximation 

In general, the computation should proceed through the converging-diverging section to the nozzle 
exit. Here, the flow is supersonic, and all of the characteristics are outrunning so that no physical condi- 
tions need be specified. However, it is often sufficient to employ a subsonic outflow boundary condition 
at the nozzle entrance by making the well-known short-nozzle approximation. 10 ’ 11 That is, the flow is 
assumed to accelerate to sonic velocity in zero length, which is equivalent to fixing the Mach number at 
the subsonic outflow boundary, M' = 0, for the single in-running characteristic. This approximation is 
strictly valid for pure longitudinal-mode oscillations only and, at worst, underestimates the nozzle damp- 
ing effect. To obtain a mean flow solution prior to making a combustor stability calculation, it is necessary 
to specify the static pressure at the nozzle entrance such that it matches the mean chamber pressure for the 
rocket engine. 
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The short-nozzle approximation can also be expressed in terms of an admittance relationship in the 


form 




e 


(41) 


where a sn is the short-nozzle admittance coefficient. It is of interest to derive an expression for this admit- 
tance coefficient using small disturbance linearization theory. 


First, note that equation (41) is equivalent to the differential form 

dp _ pdu + udp _ du i dp ( A 0 ^ 

a sn — — I • ( 4_ ) 

p pu up 

Then, consider the Mach number definition, M = u / c, the isentropic relationship for an ideal gas, 
p I p^ -k, and the acoustic velocity relationship for an ideal gas, c = sjyp / p . Taking the logarithm of 
these expressions and putting them in differential form yields 


dc du dp dp dc 1 dp dp 

c u ’ p p c 2 _ p p 


(43) 


where the short-nozzle condition, dM = 0, has been enforced. Using these relationships to elimi- 
nate du/u and dpi p in equation (42), we deduce the following value for the short-nozzle admittance 
coefficient: 


y + 1 
= ' 


(44) 
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3. VALIDATION 


As a baseline validation, it is convenient to examine the simplified case of unsteady, nonreacting, 
isentropic duct flow and to compare the numerical scheme predictions with small disturbance linear acous- 
tics theory. In this way, the fidelity of the numerical scheme can be evaluated with respect to the basic 
underlying fluid dynamic phenomena in unsteady flow. Quantitative assessment of acoustic wave driv- 
ing/damping mechanisms associated with chemical reaction processes is critically sensitive to the detailed 
physical submodeling and is ultimately dependent on empirical validation. Identification and verification 
of these driving and damping mechanisms are, of course, the central objectives of combustion instability 
research. 


3.1 Linearized Small Perturbation Acoustics 


For constant area duct flow, the continuity and momentum perturbation equations for an inviscid 
fluid may be combined to obtain an acoustic wave equation in the form 12 


(l — M 2 



i ay 

c dt 2 


2 M d 2 p' 
c dxdt 


(45) 


where the usual linearization definitions apply. Note that the cross derivative term in equation (45) arises 
as a consequence of the nonzero mean flow velocity, u ^ 0, in the duct. Considering a complex harmonic 
solution of the form p' <=c e l,LlA e li2t in equation (45), we obtain the dispersion relationship 


** ~ = o , 

’ C c 


(46) 


which is quadratic in terms of the wave number p . Thus, the characteristic wave numbers are 


Q -Q 

^ ~ (l-M)c ’ ~~ (l + M)c 


(47) 


and the complex harmonic solution takes the general form 

^ = K l [e i ^ x + K 2 e^ 2X Y Qt , (48) 

where K\ and Kj are undetermined constants. 

Additional constraints follow from a consideration of the characteristic form of the Euler 
equations: 13 
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( 49 ) 


ds ds 

— + u — = 0 
dt dx 


and 


dt 


u + 


7-1 


+ {u + c ) — 
dx 


u + 


7-1 


= 0 


(50) 


d_ 

dt 


u - ■ 


7-1 


+( u ~ c )y 

dx 


u - ■ 


7-1 


0 


( 51 ) 


which express the propagation of entropy along the streamline Cq as defined by dx/dt = u and the 
propagation of pressure waves, as described by the Riemann variables W 2 = u + 2c/(y - 1) and 
W 3 = u — 2 c/(y - 1), along the characteristics C + and C_, as defined by dx/dt = u + c and dx/dt - u- c, 
respectively. For instance, linearization of the entropy conservation relationship about the mean value 
implies s' = 0 (in the absence of flow discontinuities), and linearization of the Riemann characteristic 
equations yields 


and 


dt 


u + ■ 


7-1 


+ [(« +w') + (c +C 7 )] — 
ox 


u + ■ 


7-1 


= 0 


d_ 

dt 


u - ■ 


7-1 


+ [(« +u')-(c + c / )l ~- 
dx 


u - ■ 


7-1 


0 , 


from which complex harmonic solutions may be deduced: 


and 



yM \ 


K 2 e^ Ae iQt 


(52) 


(53) 


(54) 


C ^_ = K \{Y 0 / g ifi i x + K J^_ X ) jOt . ( 55 ) 

c 2 y \ ’ 

In order to facilitate further development of the generalized solutions, it is useful to introduce the 
acoustic impedance function z as a parameter using the standard definition 



( 56 ) 
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Application of this impedance condition at the duct inlet implies ( p' j p}f (u'/u )| = Zi, P and substitution 

of the complex harmonic solutions yields the result 


_ yM Zj n 
yM-Zin 


(57) 


The impedance condition may also be applied at the duct exit to obtain [p' / p)/[u'/u 
substitution of the complex harmonic solutions gives the additional relation, 



= where 


c *(m?.-Mi) _ yM z ex 1 ( 5 g) 

yM-z ex K 2 

Thus, equations (57) and (58) may be combined to obtain a relation between the characteristic wave num- 
bers and the inlet/outlet acoustical impedances 


Mi) h _ 

" yM + z ex "j 

f yM-z it A 


v yM — z ex j 

v yM + Zi n j 


(59) 


Taking the natural logarithm of this equation yields the more convenient form 


fe-Ah) ld= ln 


n = 




2L 


/In 


(yM + z ex "j 


_v yM ~ z ex J 

v yM + Zi n J ~ 

, it is possible to develop a 

( y M + Zex^ 


_v y^ — z ex j 

l yM + z in ) 


- i2nn , 


(60) 


+ 2 nn 


(61) 


Now, recall that the complex harmonic solution for acoustical disturbances at any fixed location in 
the duct is proportional to e 1 ^ 1 = e icot , where an amplification coefficient, X, and oscillation frequency, 
co, have been defined such that 


c M - 1 
i£2 = X + ico = — -in 

2 l d 


f yM + z ex "j 

(yM-z ln y 

_v yM — z ex y 

v yM + Zi n j _ 


+ l - 


1 -M 2 ) 


Id 


-nn 


(62) 


Thus, the amplification coefficient for acoustical disturbances in the duct is defined by the generalized 
expression 


X = 


c M 2 



2 l d 


ln 

( yM + Z ex ' 

( yM - Zi„ 


v yM — z ex y 

v yM + <-in J 


(63) 
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and the period of oscillation, P, for a given acoustical mode, n, is given by 


P = -r=—= • ( 64 ) 

f a> nc\\-M l \ 

Note that the value of the amplification coefficient governing acoustical growth/decay rates is 
dependent on the impedance functions for the duct inlet and exit. These parameters are determined from 
the application of appropriate boundary conditions, as described below. 

For subsonic inflow, there is only one outgoing characteristic and two physical boundary condi- 
tions must be specified. Thus, we impose constant mass flux and constant entropy conditions at the duct 
inlet: 


d M Lo =0 => d[ln(pu)] | JC=0 =rf[ 1 n(p) + ln(«)] | x=0 = 0 (65) 


and 


d 



= 0 

x=0 


from which we deduce the following 
forms: 


In 


(p/p 7 ) =d[ln(p)-yln(p)] 


U=o 


0 , 


( 66 ) 


differential relationships and their equivalent linearized acoustic 


and 


du _ dp u' p' 
u pup 


dp dp 

— = r— 

p p 



(67) 


( 68 ) 


Note that equation (68) applies not only at the inlet but throughout the duct since entropy must remain 
invariant along the entire flow path as previously deduced as a consequence of equations (49)-(51). Com- 
bining equations (67) and (68) yields an expression defining the impedance function at the duct inlet, 


P_ 

P 


x=0 


u 

-y— 

u 


Zin I 


x=0 


(69) 


For subsonic outflow, there are two outgoing characteristics and it is only necessary to specify one 
physical boundary condition. In this case, we introduce an acoustical admittance condition for mass flux 
at the duct exit: 


pit 

pu 


x=L 


x=L 


(70) 
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from which we deduce the following differential relationship and its equivalent linearized acoustic form: 


du dp dp u p p 

— + — = a— => — + ^ = a^r ■ 
up p u p p 


(71) 


Eliminating p' / p using equation (68), which applies throughout the duct, yields an expression defining 
the impedance function at the duct exit: 


P_ 

P 


y u 


x=L 


1 -ya u 


Zex 


r 


x=L 


1 -ya 


(72) 


The working relation for the acoustic amplification coefficient follows from elimination of the 
inflow and outflow impedances in equation (63) using the results of equations (69) and (72). 


3.2 Evaluation of Numerical Methodology 


The numerical methodology is evaluated and validated by comparing CFD predictions with small 
disturbance linear acoustics theory for unsteady, isentropic duct flow. Here, a small amplitude first-mode 
pressure perturbation is imposed on the computed mean duct flow, and the computation is restarted to 
observe the temporal evolution of this acoustical disturbance for a given outflow admittance coefficient. 
By arbitrarily specifying the exit-plane admittance coefficient, it is possible to obtain growing or decaying 
acoustic waves as desired. The resulting stability attributes, as defined by the amplification coefficient and 
period of oscillation, may then be directly compared with linear acoustics theory. For simplicity, the 
physical properties of air were utilized for all calculations. The physical conditions for the baseline duct 
acoustics problem were = 0.5 m, M = 0.2, p = 1 atm, and T = 300 K. 

A series of computations were first carried out to examine the influence of grid density on numer- 
ical performance. These calculations were based on a fixed exit plane admittance value, a =-0.1, for 
which small-disturbance linear acoustics theory predicts A=19.56 s -1 and P = 3.00 ms (i.e., / = 333 Hz) 
for the fundamental mode. The relatively large value for A implies rapidly increasing acoustic wave 
strength culminating in the development of large amplitude acoustic shocks. The numerical calculations 
were performed with 50, 100, 150, and 200 grid points using El and CN time integration with the convec- 
tive fluxes evaluated using 2CD, 2UP, and 3UP schemes. The temporal evolution of pressure at the 
outflow boundary was then used to compute the amplification coefficient and oscillation period during 
the early growth period when the acoustic amplitudes were small. The results are summarized in 
tables 1 and 2. 


Inspection of these tabulations reveals that the grid density has a strong effect on the predicted 
amplification coefficient but an extremely weak effect on the oscillation period. It is also evident that the 
first-order El technique yields poor stability predictions, irrespective of flux differencing scheme, and 
even results in wave decay for grossly coarse grids. The second-order CN scheme, on the other hand, per- 
forms much better. Although it slightly under predicts wave growth rate, the method is probably adequate 
for determining approximate stability limits. The best all-around results were obtained when using CN 
time integration with an upwind flux differencing scheme. 
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Table 1. Computed amplification coefficient, A (1/s). (Linear theory: a=-0. 1 , 


^linear = 19 - 56 s and linear = 3 mS ‘) 

^grid 

EI/2CD 

EI/2UP 

EI/3UP 

CN/2CD 

CN/2UP 

CN/3UP 

50 

-38 

-40 

-40 

12.5 

10 

10 

100 

-7.1 

-9.2 

-9.1 

14.3 

15.1 

15.1 

150 

1.8 

1.3 

1.3 

14.8 

15.9 

15.9 

200 

6.3 

6.1 

6 

15.6 

16.4 

16.4 

Table 2. Computed oscillation period, P (ms). (Linear theory: a = - 

^linear = 19 - 56 ^ and ^linear = 3 ms 0 

-0.1, 

N 

'"grid 

EI/2CD 

EI/2UP 

EI/3UP 

CN/2CD 

CN/2UP 

CN/3UP 

50 

2.98 

2.98 

2.98 

2.92 

2.98 

2.98 

100 

2.99 

2.99 

2.99 

2.98 

2.99 

2.99 

150 

2.99 

2.99 

2.99 

2.99 

2.99 

2.99 

200 

2.99 

2.99 

2.99 

2.99 

2.99 

2.99 


Comparison of small disturbance linear acoustic theory with the CN/upwind differencing CFD 
method indicates the need for a minimum grid density to obtain reliable predictions. In general, this can 
be expressed as a minimum number of grid points per wavelength for the fundamental longitudinal 
mode: 




Ngrid 

2 Id 


>100 nodes/wavelength . 


(73) 


As an example, figure 3 shows the predicted pressure fluctuations at the duct exit for a 200-grid-point 
CFD calculation superimposed with the amplitude envelope from small disturbance linear acoustic theory. 
To better illustrate the comparison, figure 4 presents the same data over a smaller time interval with the 
corresponding power spectrum inset. Here, the CFD results closely follow the classical exponential growth 
rate predicted by linear theory until the oscillations become large enough to induce significant nonlinear 
effects. Beyond this point, the wave front gradually steepens into an acoustic shock as the amplitude 
growth rate declines and a limit cycle oscillation is established. This example also serves to illustrate the 
nonoscillatory, high-resolution features of the numerical scheme. 


As further validation, it is of interest to examine the CFD predicted amplification coefficient as a 
function of the outflow admittance coefficient in comparison to linear acoustic theory. These results are 
shown in figure 5 over a relevant range of admittance values. In general, the CFD methodology slightly 
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under predicts the wave amplification rate over the full range of admittance values. This behavior is more 
pronounced when the admittance magnitudes become large and nonlinear mechanisms have an apprecia- 
ble effect. Because the CFD predicted instability growth rates err on the conservative side, it is believed 
that the method is sufficiently accurate for studying alternative combustion models and their resulting 
impact on system stability. 
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4. APPLICATION 


Some sample combustor stability calculations are now examined for illustrative purposes. The 
representative results are intended to provide concrete utilization examples and demonstrate how combus- 
tion models may be used to study instability characteristics. The first step in evaluating combustor stabil- 
ity characteristics is to compute a mean flow solution using the two-zone combustion model with all 
sensitivity indices set to zero. This also requires specification of the mean chamber pressure at the nozzle 
entrance. Once a fully converged steady-state solution has been obtained, the calculation may be restarted 
using nonzero combustion response indices. In this case, the short-nozzle admittance approximation is 
invoked for the outflow boundary condition. From a practical perspective, rocket combustor instabilities 
arise from broadband disturbances in the chamber, and each mode will grow or decay according to its own 
response attributes. Here, the initial disturbance is simply provided by a pressure perturbation of the first 
longitudinal-mode oscillation. 


4.1 Generic Rocket Chamber 

For the example calculations, representative physical conditions and propellant properties were 
specified for a generic liquid-propellant rocket engine. The geometric dimensions for this generic rocket 

O 

chamber were l c = 0.5 m and A in = A ex = 0.05 m , and the chamber conditions were p ex = 34 atm, 
Tflame = 2, 800 K, and m tot = 50 kg/s. The thermodynamic properties of the burned propellant were taken 
to be y = 1.2 and C p = 1, 800 rn 2 /s“ ■ K. All calculations were performed on a 100-grid-point mesh using 
the CN/2UP algorithm. 

Two distinct cases were examined using the two-zone combustion model. Case I assumed that the 
entire combustion process was collapsed on the injector face such that all propellant entered the chamber 
in gaseous form at a mean temperature corresponding to the specified flame temperature, 
T s = 7j lame = 2, 800 K. The mean transpiration Mach number at the inflow boundary was therefore 
adjusted to obtain the proper total mass flow rate of the rocket (i.e., M s = 0.243). Case II assumed one- 
fifth of the propellant was reacted in a collapsed combustion zone at the injector face while the remainder 
was injected as a well-dispersed spray with zero atomization length. In this case, the mean transpiration 
temperature and Mach number were taken to be T s = 1,000 K and M s =0.03, respectively. The spray 
droplets had an SMD of 100 pm, a density of 750 kg/m 3 , and entered the chamber with a velocity of 
75 m/s. Upon injection, the physical droplets were grouped into numerical parcels represented by a uni- 
form spatial distribution with a half- width equal to the computational grid spacing. The temporal behavior 
of the entire collection of physical droplets was then determined from the history of a single droplet as it 
was transported and reacted according to the distributed combustion model, assuming a mean droplet burn 
rate of 1 mm-/s. 
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4.2 Stability Limits 


4.2.1 Case I— Fully Collapsed Combustion 

Because the entire combustion process is fully collapsed onto the injector face for case I condi- 
tions, the burned gas simply enters the chamber at the specified flame temperature with a velocity that 
satisfies mass conservation. The absence of distributed combustion leads to a trivial mean flow solution 
through the chamber. The stability characteristics of the combustor were then examined by invoking the 
short-nozzle admittance approximation and restarting the calculation from the mean flow solution with a 
1 -percent initial pressure disturbance to the first longitudinal-mode oscillation. The mass transpiration rate 
of the collapsed combustion zone was made sensitive to pressure fluctuations by imposing nontrivial com- 
bustion response attributes using the two-parameter phenomenological model defined in equation (30). 

For demonstration purposes, we arbitrarily set (3 S = 0 and varried a s to determine system stability 
limits. When a s = 0.9, for instance, the combustor exhibited a stable response to the imposed disturbance, 
as illustrated in the decaying pressure waveform of figure 6. However, if the combustion response index 
was slightly increased, such that a s = 1, the combustor became unstable and the pressure wave amplitude 
was found to grow, as shown in figure 7. As the combustion response index became even larger, the ampli- 
fication rate continued to increase as illustrated for a s = 1.1 and a s =1.2 in Figures 8 and 9, respectively. 
The resulting Case I combustor stability characteristics are summarize in figure 10, which shows the pre- 
dicted amplification coefficient as a function of the pressure-sensitive combustion response index, a s . 
The stability limit for this particular case is a s ~ 0.95. 

4.2.2 Case II— Partially Collapsed Combustion 

For case II conditions, the combustion process is partially collapsed onto the injector face with 
20 percent of the total mass flow entering the chamber through a transpiring reactive boundary. The 
remaining 80 percent of the propellant is injected as a spray and consumed according to the distributed 
combustion model. The converged mean mass flow and Mach number distributions for case II conditions 
are shown in figure 11. Again, stability characteristics of the combustor were examined by invoking the 
short-nozzle admittance approximation and restarting the calculation from the mean flow solution with a 
1 -percent initial pressure disturbance to the first longitudinal-mode oscillation. The difference from case 
I being that the mass transpiration rate of the collapsed combustion zone was made insensitive to pressure 
fluctuations (i.e., a s = (5 S = 0) whereas the distributed combustion process was made sensitive to pressure 
and velocity fluctuations by imposing nontrivial combustion response attributes using the two-parameter 
phenomenological model defined in equation (40). 

For demonstration purposes, we arbitrarily set f5 v = 0 and varried a v to determine system stability 
limits, assuming zero phase lags {^ a = T fi = 0)- Figure 12 shows the unstable response of the combustor 
to the initial disturbance when a v = 1. Here, combustion pressure coupling was strong enough to over- 
come damping effects and acoustic nozzle losses to drive instability. The pressure wave amplitude was 
observed to grow exponentially until nonlinear effects lead to the formation of a limit cycle oscillation 
with a peak fluctuation >10 percent of the mean pressure. The resulting limit cycle oscillation is periodic 
but not perfectly sinusoidal with nonlinearities clearly evident in the waveform. This observation is rein- 
forced by the spatial pressure profiles within the unstable combustor, as shown in figure 13 at various time 
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intervals, where the pressure wave is impinging the injector face and being reflected back toward the 
nozzle. The resulting case II combustor stability characteristics are summarize in figure 14, which shows 
both the predicted amplification coefficient and the limit cycle peak-to-peak pressure fluctuation as a 
function of the pressure-sensitive combustion response index, a v . The stability limit for this particular 
case is a v ~ 0.7. 


22 


5. CONCLUSIONS 


A computational method was developed for the analysis of longitudinal-mode liquid-propellant 
rocket combustion instability based on the unsteady, quasi-one-dimensional Euler equations with appro- 
priate source terms. These equations were integrated in time using second-order, high-resolution, charac- 
teristic-based, flux differencing spatial discretization with Roe averaging of the Jacobian matrix. A 
two-zone combustion model was introduced where the propellant is either reacted in a collapsed combus- 
tion at the injector face or in a distributed combustion process with interphase transport derived from a 
Lagrangian treatment of representative spray droplets. The local combustion rate was made sensitive to 
pressure and velocity fluctuations through the simple introduction of proportional combustion response 
indices. It was also necessary to properly implement reflectance/admittance boundary conditions for 
impinging acoustic waves in such a way that damping mechanisms would be accurately represented. 

For baseline validation, it was convenient to compare CFD predicted wave amplification rate and 
oscillation frequency with small disturbance linear acoustics theory for unsteady, nonreacting, isentropic 
duct flow. The major findings of that validation effort were as follows: (1) The best all around CFD per- 
formance was obtained when using CN time integration with an upwind flux differencing scheme and 
(2) reliable CFD predictions could only be obtained when the computational mesh exceeded a minimum 
grid density per wavelength, N \ = 100. Although the second-order CN scheme under predicts wave growth 
rate to a slight degree, the method is generally adequate for determining system stability limits based on 
available combustion process models. 

Although the current development has been confined to an over-simplified linear combustion 
response model, more comprehensive physical submodeling can be readily implemented as desired by the 
analyst. In fact, it appears that the most important use of CFD tools for rocket combustion instability is to 
serve as a research test-bed for investigating the effect of alternative physical submodels on underlying 
processes. For illustrative purposes, some sample stability calculations were carried out for a generic com- 
bustor configuration. The objectives of this basic exercise were to demonstrate computer code utilization 
procedures for instability calculations and investigate the effect of our simplified collapsed and distributed 
combustion process models on instability characteristics. Using this approach, it was shown how the com- 
putational methodology could be used to directly determine linear stability limits as well as reveal impor- 
tant nonlinear effects— particularly as they relate to the development of steep fronted acoustic waves and 
long-term limit cycle oscillations. 

In conjunction with experimental data, much could be learned about the chemico-physical nature 
of combustion driven instabilities through CFD analyses. With time, as modeling becomes more refined 
through interplay between experiment, theory, and computation, it should be possible to evolve a predic- 
tive capability that would directly aid and support engine design and development activities. It is hoped 
that the computational framework developed herein will serve as a meaningful contribution towards that 
goal. Beyond submodel refinements and essential validation efforts, the next logical evolutionary step 
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would be to extend the framework to multiple dimensions for inclusion of tangential-mode instabilities, 
which are known to play a central role in the most dangerous and destructive forms of liquid-propellant 
rocket resonant combustion. 
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Figure 1. Generalized schematic of a conventional liquid-propellant rocket engine. Liquid 

propellant enters the chamber through an injector face plate at the head of the engine, 
and the combusted gases are expanded through a converging-diverging nozzle. 
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Figure 2. Schematic of computational domain indicating atomization region, collapsed 
and distributed combustion model zones, and nozzle entrance plane. 
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Figure 3. CFD-predicted exit-plane pressure oscillations superimposed with the amplitude 
envelope from small disturbance linear acoustics theory. Outflow admittance 
coefficient is a = -0. 1 . 
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Figure 4. Exploded view of CFD-predicted pressure oscillations superimposed with the amplitude 
envelope from small-disturbance linear acoustics theory with the corresponding power 
spectrum inset. Outflow admittance coefficient is a = -0. 1 . 


Duct Acoustics (air) 
L = 0.5 m 


Amplification Coefficient (1/s) 



30 


Figure 5. CFD-predicted amplification coefficient as a function of the outflow admittance 
coefficient in comparison to linear acoustic theory. 


Collapsed Combustion Model 



31 


Figure 6. Case I combustor response to 1 percent pressure disturbance with a s = 0.9. 
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Figure 7. Case I combustor response to 1 percent pressure disturbance with a 
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Figure 8. Case I combustor response to 1 percent pressure disturbance with a s = 1.1 . 
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Figure 9. Case I combustor response to 1 percent pressure disturbance with a s = 1.2. 
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Figure 10. Predicted amplification coefficient as a function of the pressure-sensitive 
response index for case I conditions. 
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Figure 11. Mean mass flow and Mach number distributions in the combustion chamber 
for case II conditions. 
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Figure 12. Case II combustor response to 1 percent pressure disturbance with a 
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shown at various time intervals. 
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Figure 14. Predicted amplification coefficient and limit cycle peak-to-peak pressure fluctuation 
as a function of the pressure-sensitive response index for case II conditions. 
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